#include "cem/material_law.h"
#include "cem/electromagnetics_system.h"
#include "cem/constants.h"
#include "cem/material_utils.h"

#include <array>
#include <cmath>

namespace cem {

using Tensor3x3 = std::array<std::array<double, 3>, 3>;

static inline Tensor3x3 make_isotropic(double v)
{
    return Tensor3x3{{{v, 0.0, 0.0}, {0.0, v, 0.0}, {0.0, 0.0, v}}};
}

std::array<std::array<double, 3>, 3>
LegacyMapMaterialLaw::get_sigma_tensor(int subdomain_id, const ElectromagneticsSystem & system) const
{
    Tensor3x3 sigma{{{0,0,0},{0,0,0},{0,0,0}}};
    const auto & aniso_cond = system.material_manager().area_aniso_conductivity();
    if (auto it = aniso_cond.find(subdomain_id); it != aniso_cond.end())
    {
        const auto & s = it->second;
        sigma = cem::material::build_tensor_from_principal(
            s.principal_1, s.principal_2, s.principal_3,
            s.strike_deg, s.dip_deg, s.slant_deg);
    }
    else
    {
        const auto & cond_map = system.material_manager().area_cond_map();
        const double cond_scalar = (cond_map.count(subdomain_id) ? cond_map.at(subdomain_id) : 0.0);
        sigma = make_isotropic(cond_scalar);
    }
    return sigma;
}

std::array<std::array<double, 3>, 3>
LegacyMapMaterialLaw::get_epsl_tensor(int subdomain_id, const ElectromagneticsSystem & system) const
{
    Tensor3x3 epsl{{{0,0,0},{0,0,0},{0,0,0}}};
    const auto & aniso_eps = system.material_manager().area_aniso_dielectric();
    if (auto it = aniso_eps.find(subdomain_id); it != aniso_eps.end())
    {
        const auto & s = it->second;
        epsl = cem::material::build_tensor_from_principal(
            (s.principal_1) * cem::constants::vacuum_permittivity,
            (s.principal_2) * cem::constants::vacuum_permittivity,
            (s.principal_3) * cem::constants::vacuum_permittivity,
            s.strike_deg, s.dip_deg, s.slant_deg);
    }
    else
    {
        const auto & epsl_map = system.material_manager().area_epsl_map();
        const double epsl_scalar = (epsl_map.count(subdomain_id) ? epsl_map.at(subdomain_id) : 1.0);
        const double epsl_eff = (epsl_scalar) * cem::constants::vacuum_permittivity;
        epsl = make_isotropic(epsl_eff);
    }
    return epsl;
}

std::array<std::array<double, 3>, 3>
LegacyMapMaterialLaw::get_inv_mu_tensor(int subdomain_id, const ElectromagneticsSystem & system) const
{
    Tensor3x3 inv_mu{{{0,0,0},{0,0,0},{0,0,0}}};
    const auto & aniso_mag = system.material_manager().area_aniso_magnetic();
    if (auto it = aniso_mag.find(subdomain_id); it != aniso_mag.end())
    {
        const auto & s = it->second;
        auto inv_mu_r = cem::material::build_inv_tensor_from_principal(
            s.principal_1, s.principal_2, s.principal_3,
            s.strike_deg, s.dip_deg, s.slant_deg);
        const double inv_mu0 = 1.0 / cem::constants::vacuum_permeability;
        for (int i = 0; i < 3; ++i)
            for (int j = 0; j < 3; ++j)
                inv_mu[i][j] = inv_mu_r[i][j] * inv_mu0;
    }
    else
    {
        const auto & mu_map = system.material_manager().area_mu_map();
        const double mu_r_eff = (mu_map.count(subdomain_id) ? mu_map.at(subdomain_id) : 1.0);
        const double mu_r_safe = (std::fabs(mu_r_eff) > 0.0 ? mu_r_eff : 1.0);
        const double inv_mu_scalar = 1.0 / (cem::constants::vacuum_permeability * mu_r_safe);
        inv_mu = make_isotropic(inv_mu_scalar);
    }
    return inv_mu;
}

} // namespace cem